% Scripts that read res files, assemble result matrix ifft the date back to
% time series
%simulation for Chile tsunami em signals
% Hx1  = -(Htr+j*Hti)*100*4.*pi;
% Hy1  =  (Hpr+j*Hpi)*100*4.*pi;
% Hr1  =  (Hrr+j*Hri)*100*4.*pi;

%%
Hr = zeros([361,100,396]);
prd = 1./(fftfrq(720,1/(60/(24*3600))));
prd = [4000. prd]; % The mean flow period -> 4000 days (it should have been DC but 4000 days OK)!
base_dir = 'D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\resfiles\';
for i = 1 : 200,
display(['    ' sprintf('%10.5f',prd(i)) '     ' sprintf('tsu%d',i) ' tsunami25.model ' sprintf('tsu%d',i) '.source' ' 1']);
[Hx1,Hy1,Hr1] = read_chile_res([base_dir sprintf('tsu%d',i) '.res']);
Hr(i,:,:) = Hr1(:,805:1200);
end;
%Note. Compared ts synthesized from 131 harmonics to that with 200. No
% noticeable difference. hence no need of going to harmonics higher than
% 131 the saved data files are all from the first 131 harmonics
%

%%

%Hr(132:361,:,:) = 0; %filling up the higher harmonics for which the no simulations were done.
% Hp(241:301,:,:) = 0;% this is justified, since the higher harmonics of t flow are
% weaker.
Hr(201:361,:,:) = 0;

for i = 1:100,
    for j = 1:396,
         Hr_t(:,i,j) = real(ifft([squeeze(Hr(:,i,j)).' fliplr(squeeze(Hr(2:end-1,i,j))')])); % Hr(2:end-1,i,j) because the DC to be omitted and the last value is the same (fft symm)
%          Hp_t(:,i,j) = real(ifft([squeeze(Hp(:,i,j)).' fliplr(squeeze(Hp(2:end-1,i,j))')]));
%          Ht_t(:,i,j) = real(ifft([squeeze(Ht(:,i,j)).' fliplr(squeeze(Ht(2:end-1,i,j))')]));
    end;
end;

%%load
fday_simulation=datenum(2010,2,27,06,34,14+time_data);
plot(fday_simulation,squeeze(Hr_t(:,50,204)));
datetick
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ipc_data\ipc20100227vsec.sec.mat;
ipc_z = data.DATA(:,3);
fday = data.FDAY;
b1 = min(fday):(1/24):max(fday);
sp=spline(b1,ipc_z(1:100:end)'/spline(b1,eye(length(b1)),fday(1:100:end)'));
ipc_z_f=ipc_z-ppval(fday,sp);
ipc_h = data.DATA(:,1);
sp=spline(b1,ipc_h(1:100:end)'/spline(b1,eye(length(b1)),fday(1:100:end)'));
ipc_h_f=ipc_h-ppval(fday,sp);
plot(fday,ipc_h_f)

%% Test
%% 
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\results_matfiles\Hr_spectra_time_series Hr_t;
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\modeling\results_matfiles\Easter_Island_ts;
%%
% Easter Islands -109.41 and -27.17
% Grid points 50,203

for i = 50:50,
    for j = 203:203,
        plot(fday,ipc_z_f,'r');
        hold on;
        plot(fday_simulation,squeeze(Hr_t(:,i,j)),'k','LineWidth' , 2);
        title(sprintf(' i = %d j = %d', i,j));
        axis([min(fday_simulation) max(fday_simulation)  -1 1]);
        pause;
        hold off;
    end;
end;